---
title: 'Both: Cox model'
author: 'David Hume'
date: '`r format(Sys.time(), "%d %B %Y")`'
output: 
   html_document:
    toc: true
    toc_depth: 3 
editor_options: 
  chunk_output_type: console
---
---

```{r setup, message=FALSE, include = FALSE}
library(tidyverse)
library(data.table)
library(Hmisc)
library(fasttime)
library(lubridate)
library(stringr)
library(cowplot)
library(car)
library(stargazer)
library(lfe)
library(plm)
library(knitr)
library(pryr)
library(DescTools)
library(broom)
library(survival)
library(survminer)
options(width=200, scipen=10)
knitr::opts_chunk$set(echo = TRUE, cache = FALSE, warning = FALSE, 
                      message = FALSE, cache.lazy = FALSE)

```

<div style="line-height: 2em;">
<font size="4"> 

---

```{r echo = FALSE, eval=TRUE}

Time <- Sys.time()

# Set-up

# Switches & setting of variables

# Data set either "All8pc", " "All1pc", "Test", "Render"
analysis.1 <- "Render"
# Filetype for output
# "latex" 
# "html" - inserting into html [htmltools::includeHTML()] or Word docs [Insert/[Text]/Object/Text from file]
# "text" - viewing on screen
filetype.1 <- "Render"


analysis <- ifelse(analysis.1 =="Render", readRDS("./Data/Analysis.rds"), analysis.1)
filetype <- ifelse(filetype.1 =="Render", readRDS("./Data/Filetype.rds"), filetype.1)


if (filetype=="latex") {
  out <- "tex"

  } else if (filetype=="html") {
      out <- "html"

    } else if (filetype=="text") {
      out <- "csv"

      } else {
        stop("latex, html or text")
}



# Set variables
fig = 0 # Counter for figures

```

```{r cache=FALSE, eval=TRUE, echo=FALSE}

Keep_cols <- c("TUI", "Property_Id", "Type", "New", "Quality100000", "Region", "date.val", "sold.dummy", "IdxPriceA", "PriceA",  "Peak_price")




#if (analysis.1!="Render") {
  if (analysis=="All") {
    f <- function(x, pos) subset(x, Type == Type) # i.e. reads in all data because condition always TRUE
      Resold.indexed <-  read_csv_chunked("./Data/Resold.indexed.All.csv", DataFrameCallback$new(f), chunk_size = 2000000)
      Resold.indexed <- as.data.table(Resold.indexed)
      Resold.indexed.regress <- Resold.indexed[, ..Keep_cols]
      Resold.indexed <- NULL
      
      } else {
        Resold.indexed.regress <-fread(paste0("./Data/Resold.indexed.", analysis, ".csv"), select = Keep_cols)
  }
#}



Resold.indexed.regress <- Resold.indexed.regress[,date.val:=as.Date(date.val, origin ="1970-01-01")]



cat("DATA:", analysis)

```





```{r cache=FALSE, eval=FALSE, echo=FALSE}

# Zhang et al (2018) Time-varying covariates and coefficients in Cox regression models

library(survsim) # Need to update version of R, works ok on laptop
library(survival)

N=100 #number of patients
set.seed(123)
df.tf<-simple.surv.sim(#baseline time fixed
  n=N, foltime=500, dist.ev=c('llogistic'), anc.ev=c(0.68), beta0.ev=c(5.8), anc.cens=1.2,
  beta0.cens=7.4, z=list(c("unif", 0.8, 1.2)), beta=list(c(-0.4),c(0)), x=list(c("bern", 0.5), c("normal", 70, 13)))
names(df.tf)[c(1,6,7)]<-c("id","grp","age")



set.seed(123)
nft<-sample(1:10, N,replace=T)#number of follow up time points crp<-round(abs(rnorm(sum(nft)+N,
crp<-round(abs(rnorm(sum(nft)+N, mean=100,sd=40)),1) 
time<-NA
id<-NA
i=0
for(n in nft){ 
  i=i+1
  time.n<-sample(1:500,n) 
  time.n<-c(0,sort(time.n)) 
  time<-c(time,time.n) 
  id.n<-rep(i,n+1) 
  id<-c(id,id.n)
}
df.td <- cbind(data.frame(id,time)[-1,],crp)


df<-tmerge(df.tf,df.tf,id=id, endpt=event(stop,status))
head(round(df))


df <- tmerge(df,df.td,id=id, crp=tdc(time,crp))
head(round(df),10)

fit.tdc <- coxph(Surv(tstart,tstop,endpt)~ grp+age+crp+cluster(id),df)
fit.tdc



```




```{r cache=FALSE, echo=FALSE, results = FALSE}

# Equivalent code incuded within regression variables function
Resold.indexed.regress<-Resold.indexed.regress[!is.na(Peak_price),]
Resold.indexed.regress<-Resold.indexed.regress[, F_unrealgain :=round(IdxPriceA-Peak_price) # Calculates unrealised gain or loss - price paid vs index
                                       ][, `:=`(F_unrealgainAdj = Winsorize(F_unrealgain/100000, probs = c(0.01, 0.99), na.rm = TRUE),
                                                F_unrealgainYN = F_unrealgain >=0,
                                                Price_unrealgainAdj = Winsorize((IdxPriceA-PriceA)/100000, probs = c(0.01, 0.99), na.rm = TRUE),
                                                Price_unrealgainYN = ifelse(IdxPriceA-PriceA>=0, TRUE, FALSE))
                                         ][, `:=`(F_unrealgain.pos = ifelse(F_unrealgainYN,F_unrealgainAdj,0),
                                                  F_unrealgain.neg = ifelse(!F_unrealgainYN,F_unrealgainAdj,0),
                                                  Price_unrealgain.pos = ifelse(Price_unrealgainYN,Price_unrealgainAdj,0),
                                                  Price_unrealgain.neg = ifelse(!Price_unrealgainYN,Price_unrealgainAdj,0))
                                           ]

Resold.indexed.regress <- Resold.indexed.regress[order(Property_Id, date.val)][,tstart:=rowid(TUI)][,tstop:=tstart+1]






```


```{r}
head(Resold.indexed.regress[,.(Property_Id, TUI, date.val, tstart, tstop, sold.dummy)], 10)

fit.coxph1 <- coxph(formula = Surv(tstart, tstop, sold.dummy) ~ Price_unrealgainYN + strata(Property_Id) + cluster(Property_Id), data = Resold.indexed.regress)

summary(fit.coxph1)


fit.coxph2 <- coxph(formula = Surv(tstart, tstop, sold.dummy) ~ F_unrealgainYN + strata(Property_Id) + cluster(Property_Id), data = Resold.indexed.regress)

summary(fit.coxph2)

fit.coxph3 <- coxph(formula = Surv(tstart, tstop, sold.dummy) ~ Price_unrealgainYN + F_unrealgainYN + strata(Property_Id) + cluster(Property_Id), data = Resold.indexed.regress)

summary(fit.coxph3)




table <- stargazer(fit.coxph1, fit.coxph2, fit.coxph3, no.space=TRUE, align=TRUE, digits = 4, digits.extra = 0,  dep.var.caption  = "",
            dep.var.labels = "Sale", object.names=FALSE, model.numbers=FALSE,  model.names=FALSE, column.labels = c("(1)", "(2)", "(3)"), 
              covariate.labels=c("Gain Since Purchase = 1", "Gain Since Peak = 1"),
              omit.stat=c("LL","ser","f", "adj.rsq", "wald", "max.rsq", "lr", "logrank"),
              float = F, omit.table.layout = "n", column.sep.width = "-15pt", style = "aer",
            table.layout="-dc#-tas-", type = "latex")


table <- str_replace_all(table, "\\^", "")
table <- str_replace_all(table, "R\\$\\{2\\}\\$", "R$^{2}$")
write.table(table[11:(length(table)-2)], col.names = F, row.names = F, quote = FALSE, paste0("./Tables/LaTeX/Cox", analysis, "Clean.tex"))

stargazer(fit.coxph1, fit.coxph2, fit.coxph3, no.space=TRUE, align=TRUE, digits = 4, digits.extra = 0,  dep.var.caption  = "",
            dep.var.labels = "Sale", object.names=FALSE, model.numbers=FALSE,  model.names=FALSE, column.labels = c("(1)", "(2)", "(3)"), 
              covariate.labels=c("Gain Since Purchase = 1", "Gain Since Peak = 1"),
              omit.stat=c("LL","ser","f", "adj.rsq", "wald", "max.rsq", "lr", "logrank"),
              float = F, omit.table.layout = "n", column.sep.width = "-15pt", style = "aer",
              type = filetype, out=paste0("./Tables/Cox", analysis, ".", out))
  






```











```{r cache=FALSE, eval=TRUE, echo=FALSE}
Sys.time() - Time

```

